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ABSTRACT 


-/ 

The  purpose  of  this  thesis  is  to  analyze  the  verti¬ 
cal  pullout  capacity  of  a  horizontal  circular  anchor  plate 
embedded  in  soil,  using  the  plasticity  theory  as  developed 
by  Sokolovskii  (1960) . 

A  mathematical  model  for  soil-structure  interaction 
was  used  to  generate  slip  line  fields  compatible  with  the 
boundary  conditions  of  the  anchor  pullout  problem.  To  see 
if  the  solution  for  a  given  set  of  boundary  conditions  was 
physically  admissible,  computer-generated  plots  of  slip  line 
fields  were  obtained.  If  the  slip  line  field  was  permissible, 
i.e.,  no  overlapping  of  characteristic  lines.  The  pullout 
capacity  was  then  calculated.  A  conventional  Coulomb  sli¬ 
ding  wedge  graphical  solution  was  performed.  The  results 
of  the  hand  solution  closely  approximate  the  computer 
results . 

The  theoretically  determined  pullout  capacities  were 
found  to  be  less  than  experimental  results  of  others.  In 
loose  sand,  at  small  depth  to  anchor  plate  diameter  (D/B) 
ratios,  the  predicted  values  were  approximately  81  percent 
of  experimental  results,  while  in  dense  sand  they  were 
approximately  50  percent  of  experimental  results. 
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CHAPTER  1 


INTRODUCTION 

The  purpose  of  this  thesis  is  to  study  the  anchor 
pullout  problem  using  the  theory  of  plasticity  and  the 
Coulomb  failure  criteria  for  soils  in  the  analysis. 

The  Problem 

The  problem  is  to  determine  the  pullout  capacity  of 
a  circular  plate  embedded  below  a  horizontal  ground  surface 
and  centrally  loaded  by  an  upward  vertical  force.  Pullout 
capacity  is  defined  as  the  maximum  vertical  load  or  pull 
required  to  remove  the  anchor  from  the  soil. 

Earth  anchors  are  used  to  provide  uplift  resistance 
for  transmission  towers,  submerged  pipelines  and  tunnels, 
mobile  homes,  aircraft  and  ocean  mooring  systems,  and  to 
develop  tieback  forces  required  to  eliminate  external  brac¬ 
ing  from  retaining  structures  and  sheeting  walls. 

Previous  anchor  pullout  determinations  are  based  on 
semi-empirical  methods  and/or  theoretical  methods  that  make 
some  limiting  assumptions.  Therefore,  there  is  a  need  for 
study  to  find  a  reliable  method  of  determining  the  anchor 
pullout  capacity. 
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Scope 

This  thesis  presents  a  theoretical,  computer  deter¬ 
mination  of  the  pullout  capacity  and  associated  slip  line 
field  for  anchors  having  various  depth  of  embedment  to  anchor 
diameter  (D/B)  ratios,  for  both  loose  and  dense  sand.  A 
description  of  the  theory  applicable  to  the  problem  is  pre¬ 
sented  first,  followed  by  a  literature  review  of  anchor 
pullout  studies,  the  solution  procedure,  a  presentation  and 
discussion  of  results,  and  finally,  conclusions  and  recom¬ 
mendations  for  future  studies. 


CHAPTER  2 


PLASTICITY  THEORY  APPLICABLE  TO  PROBLEM 

The  initial  applications  of  the  plasticity  theory 
were  mainly  in  the  field  of  metallurgy.  Calculations  based 
on  limiting  stress  fields  for  soils  were  first  reported  in 
detail  by  Sokolovskii  (1960  and  1965) .  He  used  the  method 
of  characteristics  to  derive  the  stress  fields  for  materials 
possessing  cohesion  or  friction  or  both.  It  is  Sokolovskii's 
method,  as  outlined  by  Harr  (1966)  and  Karafiath  and  Nowatzki 
(1978)  that  was  used  in  this  study.  To  understand  the 
author's  solution  procedure  better,  a  brief  description  of 
Sokolovskii's  method  is  presented  in  this  chapter. 

Fundamental  to  analyses  by  plasticity  theory  is  the 
assumption  that  the  soil  acts  as  a  rigid-perfectly  plastic 
material.  This  means  that  under  load,  the  soil  behaves  as 
illustrated  by  the  solid  line  in  Figure  1.  The  material 
does  not  deform,  regardless  of  loading  sequence,  stress 
history,  rate  of  loading,  etc.,  until  it  reaches  failure. 

This  differs  from  elastic-plastic  material  behavior  illus¬ 
trated  by  the  dashed  line  in  Figure  1.  Here  the  soil  under¬ 
goes  recoverable  deformation  before  reaching  failure.  In 
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both  cases,  deformation  at  failure  takes  place  at  a  con¬ 
stant  rate  without  change  in  stress. 
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Figure  1.  Material  behavior  (after  Karafiath 
and  Nowatzki,  1978). 

Fundamental  Concepts  and  Relationships 
The  stress  condition  at  a  point,  A,  within  a  semi¬ 
infinite  homogeneous  soil  mass  under  a  state  of  loading,  with 
horizontal  x-axis  and  vertical  z-axis,  is  shown  in  Figure  2. 
Lines  OF  and  OF'  are  the  Mohr-Coulomb  strength  envelopes. 

The  directions  of  the  lines  of  impending  rupture,  called 
slip  lines,  are  indicated  by  the  angles  y  which  are  measured 
+  and  -  from  the  direction  of  the  major  principal  stress, 
and  which  can  be  expressed  as: 

y  =  tt/4  -  4>/2  (1) 

The  orientation  of  the  major  principal  stress  with  respect 
to  the  horizontal  is  given  by  the  angle  0.  If  the  x-axis 
is  taken  along  the  horizontal,  then  0  is  the  angle  between 
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the  direction  of  the  major  principal  stress  and  the  x-axis. 
The  symbol  o  (without  subscript)  is  used  to  indicate  the 
stress  magnitude  on  the  t-o  axis  from  0,  the  point  of  inter¬ 
section  of  the  Mohr-Coulomb  envelope  and  the  o-axis,  to  the 

center  of  the  Mohr  circle,  Q.  In  general,  o  =  (o  +  o  )/ 

x  z 

2  +  i|),  where  tp  =  c  cot  4>.  Therefore,  the  stress  state  at 
any  point  can  be  expressed  as: 

ox  =  o(l  +  sin  cos  26)  -  ip 

az  =  a(l  -  sin  4>  cos  26)  -  4)  (2) 

t  =o  sin  4>  sin  26 
xz 

Figure  3  shows  the  convention  used  to  describe  the 
orientation  of  the  various  stresses  discussed  above  for  the 
point  in  question.  This  figure  demonstrates  that  the  slip 
lines  at  any  point  can  be  written  in  differential  form  as: 

gf  =  tan  (6  +  p)  (3) 

The  p(+)  sign  refers  to  " i-characteristic"  slip  lines  in¬ 
clined  at  a  slope  of  tan  (6  +  p)  from  the  x-axis  and  the 
p(-)  sign  refers  to  " j-characteristic"  slip  lines  inclined 
at  a  slope  of  tan  (6  -  p)  from  the  x-axis.  The  term  "slip 
line  field"  is  understood  to  mean  the  two  families  of  curves 
defining  the  bounds  of  regions  in  which  stress  states  cor¬ 
respond  to  those  of  incipient  plastic  failure.  The 
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directions  correspond  to  the  directions  of  shearing  stress 
causing  the  failure. 


Figure  3.  Orientations  (after  Karafiath 
and  Nowatzki,  1978). 


Assuming  that  the  x-axis  is  not  horizontal  but  is 
inclined  at  an  angle  e  to  the  direction  of  gravity,  and 
assuming  that  inertia  and  seepage  forces  are  zero,  the 
normal  and  shear  stresses  a  ,  a  ,  t  on  any  soil  element 

X  Z  X  z 

must  satisfy  the  static  equilibrium  equations: 


9t 


xz 


9z 


Y  sin  e 


9x 


xz 


9x 


Y  cos  e 


(4) 


where  y  =  effective  unit  weight. 


Substituting  Equations  2  into  Equations  4  yields  the 
basic  differential  equation  of  plastic  equilibrium  for  soil 


under  plane  stress  conditions.  Sokolovskii  modified  these 
equations  to  yield: 
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do  +  2 a  tan4>  d0  =  — ^ — - 

—  COS<t> 


[sin  ( e+(p)  dx  +  cos  ( e+<f> )  dz] 

(5) 


Numerical  Solution 

Numerical  methods  are  required  for  solution  of 
Equations  3  and  5,  for  cases  of  practical  interest.  Assuming 
that  the  x-axis  is  horizontal,  i.e.,  e  =  0,  the  equations 
become : 

For  i-characteristics  (the  slip  lines  corresponding 
to  i  =  constant) : 


dz  =  dx  tan(0  +  y)  (6a) 

do  +  2a  tan<t>  d8  =  y  (dz  +  tan  cj>dx)  (6b) 

For  j-characteristics  (the  slip  lines  corresponding 
to  j  =  constant) : 

dz  =  dx  tan ( 0  ~  y)  (6c) 

do  -  2a  tan4>  d0  =  y(dz  -  tan  <{>dx)  (6d) 

All  points  in  the  i-j  coordinate  system  of  Figure 

4a  can  be  assigned  a  position  in  terms  of  a  nodal  number 

X.  .  where  the  i  denotes  the  abscissal  reference  and  j  the 
1  •  3 

ordinate.  Associated  with  every  nodal  point  X.  .  are  the 

1  •  3 

quantities  x.  z.  ■,  o.  and  0.  ..  All  nodal  points  in 

IfJ  1 1 r  J  1 f  J 
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grid  (Figure  4a)  and  all  pertinent  information  relative  to 
each  point  can  be  put  into  tabular  form  (Figure  4b) .  The 
basic  problem  is  to  determine  the  values  of  x,  z,  o,  and  0 
at  nodal  points  of  interest  from  known  values  of  these 
parameters  at  neighboring  nodal  points. 

To  solve  Equations  6,  Sokolovskii  used  a  finite- 
difference  procedure,  and  proposed  that  they  be  replaced  by 
finite  equivalents  as  follows: 

Along  i-characteristics : 


(zi,j 

-  z  .  .  . ) 

i,  3-1 

=  (x.  .  - 

1,3 

xi,3-i’ 

tan  (0^ 

,3-1  +  u) 

(7a) 

(0i,j 

'  °i, j-l1 

+  2a .  .  . 

1,3-1 

<8i,3  - 

tan<() 

= 

Y[(zi.j  ' 

z  ,  .  , )  + 

1,3-1 

lxi,3  - 

■‘i.j-i’ 

(7b) 

tanc()] 

Along 

j-characteristics : 

(Zi,j 

‘  Zi-l,j> 

=  (x.  .  - 

if  3 

tan  (0. 

l 

-Ifj  •  U> 

(7c) 

(°i,j 

'  ai-l,j> 

- 

(ei,3  - 

tan<() 

= 

Y[(Zi,j  ' 

(xi,3  ' 

( 7d) 

tan<))] 

Equations  7  represent  four  equations  in  four  un¬ 
knowns.  Solving  these  equations  results  in  the  following 
recurrence  formulas: 
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1.1 


tan (6.  ..  . +v) -tan { 6 
1  *  l“x 


— 


(8a) 


.  =  z.  i  .  +  (x.  .  -  x.  ,  .)  tan  (0.  .  .  -  y) 

i.D  1-1*1  i*l  i-l*D  1-1*1 


(8b) 


1*1 


o .  .  ,  +  a .  ,  . 

1,3-1  1-1*1 


where 


(8c) 


C  =  (z.  .  -  z.  .  .)  -  tan<j>  (x.  .  -  x.  ..  .) 

i*l  i—l*l  i*l  i—l*l 


D  -  (Ji,j  -  zi,j-i>  *  tan*  (xi,j  - 


e.  . 

i.i 


2  tan*  (0i,j-l  +  °i-l,j> 


(8d) 


Equations  8  demonstrate  that  the  quantities  x,  z,  o, 

and  6  can  be  determined  at  a  node  X.  .if  and  only  if  they 

1,1 

are  known  at  the  two  nodes  X .  .  ,  and  X .  ,  .in  the  near 

vicinity  of  the  node.  However,  the  values  of  all  the  vari¬ 
ables  are  not  always  known  on  a  given  boundary,  nor  are 
boundaries  always  coincident  with  characteristic  lines. 
Basically,  three  types  of  boundary-value  problems  are  en¬ 
countered:  the  Cauchy  problem,  the  Goursat  problem,  and  the 
Mixed  Boundary-Value  problem. 


Cauchy  Problem 

In  this  problem  (Figures  5a  and  5b),  x,  z,  a,  and  0 
are  known  at  all  points  on  a  non-characteristic  line,  such 
as  AB,  along  which  neither  i  nor  j  is  constant.  The  problem 
is  to  determine  these  variables  along  characteristic  lines 
within  and  on  other  boundaries  of  the  triangle  ABC.  To 
accomplish  this,  AB  is  first  divided  conveniently  into  n 
sections  by  introducing  the  i  and  j  characteristics  as  shown 
in  Figure  4.  Then  the  values  of  x,  z,  a,  and  0  at  each  node 
represented  by  the  intersection  of  the  characteristics  with 
AB  ar<?  recorded  in  the  appropriate  squares  of  the  table  given 
in  Fig*  ?  5b.  These  values  occupy  the  diagonal  spaces  shown 
shaded  in  the  figure.  The  recurrence  relations  (Equations  8) 
are  used  to  determine  the  desired  values.  Proceeding  in  an 
orderly  manner  according  to  the  arrows  shown  in  the  figure, 
the  entire  region  can  be  filled  and  the  quantities  extended 
to  the  boundaries  of  the  triangle. 

Goursat  Problem 

For  this  problem  (Figures  5c  and  5d) ,  the  quantiti- 
ties  x,  z,  a,  and  0  are  known  at  all  points  on  the  charac¬ 
teristics  AB  and  AC.  As  with  the  Cauchy  problem,  the 
Goursat  problem  is  solved  by  proceeding  from  two  known 
points  to  an  unknown  point  as  shown  in  Figure  5d  until  the 
values  of  the  variables  are  known  at  all  nodes  within  the 
zone  defined  by  the  boundary  CABD. 


Figure  5.  Type  of  bounds 
1966).  —  (a,b 
problem;  (e,f) 
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Mixed  Boundary-Value  Problem 

In  this  problem,  some  of  the  quantities  x,  z,  o, 
and  0  are  known  and  some  specified  along  one  (noncharacter¬ 
istic)  boundary  and  all  of  the  quantities  are  known  along 
another  (characteristic)  boundary.  Thus  there  are  a  number 
of  variations  of  this  case.  Figures  5e  and  5f  show  one  case 
where  x,  z,  a,  and  0  are  all  known  along  the  characteristic 
AB  and  only  z  and  a  are  known  along  the  non-characteristic 
line  AC.  The  problem  is  to  determine  the  value  of  all  quan¬ 
tities  on  and  within  the  boundary  ABC.  The  recurrence 
formulas  given  by  Equations  8  will  not  directly  provide  the 
results  along  AC.  If  Equations  7c  and  7d  are  rewritten  as: 


xi,j  =  (zi,j  -  2i-l,j>  COt  (0i-l,j-^  +  xi-l,j 


(9a) 


{y[(zi  -fxi_i  i)tan4>]+(o.  .-o  .)  } 

6.  .  *  6.  ,  .+ - — L_iiJ - iiJ — iliil — - U — — 

i,D  1-1, 3  2o.  .  .  tan$ 

A  A  9  J 

(9b) 


the  resulting  expression  can  be  used  to  determine  values 
along  AC.  For  points  within  the  boundary,  where  neither  x, 
z,  a,  or  0  are  known.  Equations  8  must  be  used. 

Once  the  desired  x,  z,  a,  and  0  quantities  have  been 
obtained,  plotting  and  connecting  the  corresponding  x  and  z 
coordinates  of  the  same  characteristics  results  in  the  slip 
line  field.  The  envelope  of  these  is  the  failure  surface. 


The  corresponding  stresses  at  any  point  can  be  obtained  by 

extrapolating  the  a  and  6  values  at  the  nodal  points  and 

applying  Equations  2.  The  accuracy  of  the  results  increases 

as  the  number  of  characteristics  used  in  the  grid  increases. 

For  increased  accuracy,  Sokolovskii  recommended  that,  in 

Equations  8  and  9,  the  terms  j  be  replaced  by  ( 0 ^  ^  + 

0.  .  )/2  and  6.  .  .  by  (0.  .  .  +  0.  .  )/2.  Once  0.  .is  deter- 

mined,  0.  .  .  and  0.  .  .  are  recalculated,  as  just  described, 

i-l,D  1,3-1  J 

and  substituted  into  the  recurrence  formulas,  Equations  8 
and  9 . 


Plasticity  Approach  to 
Retaining  Wall  Problem 

The  determination  of  lateral  earth  pressures  against 
rigid  retaining  structures  is  one  of  the  most  common  problems 
encountered  in  the  evaluation  of  the  stability  of  such  struc¬ 
tures.  In  practice,  the  Coulomb  solution  is  widely  used. 

This  solution  is  a  limit  equilibrium  solution  in  which  the 
soil  is  assumed  to  fail  by  sliding  along  some  plane  inclined 
at  a  certain  angle,  uj,  from  the  backface  of  the  wall  (refer 
to  Figure  6) .  The  magnitude  of  the  angle  u  is  a  function  of 
the  internal  friction  angle,  $,  of  the  soil  behind  the  wall, 
the  wall-soil  interface  friction  angle,  6,  the  wall  backface 
inclination,  a,  with  respect  to  the  horizontal,  and  the  back¬ 
fill  inclination,  8,  with  respect  to  the  horizontal.  The 
magnitude  of  the  force  exerted  by  the  backfill  soil  on  the 
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wall  depends  upon  the  direction  of  wall  movement.  Wall 
movement  away  from  the  backfill  is  called  the  "active"  case 
and  wall  movement  into  the  backfill  the  "passive"  case. 
Conditions  imposed  by  anchor  pullout  approximate  the  passive 
case  most  closely.  In  the  passive  case  (wall  movement  into 
the  backfill)  the  soil-wall  interface  friction  angle,  6,  may 
be  positive  or  negative.  It  is  positive  for  downward  move¬ 
ment  of  the  wall  relative  to  the  backfill  and  negative  for 
upward  movement  of  the  wall  relative  to  the  backfill.  The 
(+)  and  (-)  directions  are  shown  in  Figure  6.  Therefore,  the 
wall  force,  P^,  may  act  in  any  direction. 


Figure  6.  Coulomb  solution  to  lateral  earth  pressure 

problems  (after  Nowatzki  and  Karafiath,  1981). 


Positive  6  angles  are  usually  encountered  in  conven¬ 
tional  soil  mechanics  problems  where  passive  pressures  are 
generated  by  horizontal  movement  of  the  wall  into  the  soil 
backfill,  i.e.,  the  passive  force  of  the  wall  against  the 
backfill  is  directed  downward.  In  some  cases,  however,  such 
as  soil  anchors,  bulldozer  blades  and  track  grousers,  the 
force  is  applied  at  a  negative  6  angle. 

Application  of  conventional  plasticity  theory  to  the 
problem  of  passive  earth  pressure  results  in  slip  line 
fields  that  define  the  geometry  of  the  failure  surface. 

The  failure  surfaces  shown  in  Figures  7a  and  7b  compare  well 
with  those  shown  in  literature,  as,  for  example,  by  Chen 
(1975) .  The  plasticity  solution  for  the  case  shown  in 
Figure  7c,  however,  results  in  an  overlapping  of  the  slip 
line  field.  This  overlapping  implies  that  two  different 
stress  states  exist  at  a  given  point  within  the  overlap 
region  at  the  same  time,  a  physical  impossibility.  Yet, 
plastic  zones  have  been  observed  to  form  experimentally  even 
for  the  case  <5  <  0  (Rowe  and  Peaker,  1965). 

In  plasticity  solutions  obtained  in  the  conven¬ 
tional  way,  all  three  boundary  value  problems  are  generally 
solved  in  order  to  obtain  a  complete  slip  line  field  within 
a  loaded  soil  mass.  Generally,  the  Cauchy  problem  is 
solved  to  obtain  the  passive  zone,  the  Mixed  Boundary- 
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6  >  0 
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(b) 
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assive  case  failure  zone  geometries  for  various 
alues  of  6  (after  Nowatzki  and  Karafiath,  1981) 


Value  problem  (MBV)  to  obtain  the  active  zone,  and  the 
Goursat  problem  (G)  to  obtain  the  transition  zone  between 
the  other  two. 

These  zones  and  their  related  problems  are  indicated 
in  Figure  7  for  the  case  of  passive  earth  pressure  on  a  re¬ 
taining  wall.  Zone  OAB  is  the  Cauchy  zone  and  zone  OBC  the 
Goursat  zone.  The  computation  of  x,  z,  cr,  and  0  revolves 
about  the  singular  point  0  where  o  values  are  computed  for 
specified  changes  in  0.  In  zone  ODC,  x,  z,  o,  and  0  are 
known  along  OC  after  the  Goursat  problem  has  been  solved. 

For  a  vertical  wall,  OD,  x  =  0  along  the  wall  boundary  and 
0  is  assumed.  Therefore,  this  is  a  Mixed  Boundary-Value 
problem.  Once  solved,  the  normal  pressures  along  the  wall 
can  be  determined  from  o  and  0. 

Reexamination  of  the  boundary  conditions  by  Nowatzki 
and  Karafiath  (1981)  for  this  case  resulted  in  the  following 
conclusions : 

1.  In  the  passive  case  the  wall  force  is  an  applied 
force,  the  direction  of  which  is  defined  by  the 
problem.  Therefore,  along  OD  (refer  back  to  Figure 
7)  all  four  variables,  x,  z,  o,  and  6,  are  given 
a  priori  (although  for  computational  purposes  o  and 
0  are  assumed) .  Therefore,  zone  OBD  should  be 
solved  as  a  Cauchy  problem. 
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2.  In  cohesionless  soils,  0  along  the  free  surface,  OA, 
is  undetermined  since  =  a ^  =  0.  Therefore, 

zone  OAB  should  be  solved  as  a  Mixed  Boundary-Value 
problem,  with  z  =  0  and  a  =  0  along  OA. 

Nowatzki  and  Karafiath  developed  a  computer  program 
incorporating  the  recurrence  relationships  of  Equations  8 
and  9,  to  accommodate  the  above  boundary  conditions.  There 
was  found  to  be  excellent  agreement  between  the  computed 
and  observed  failure  surfaces  of  cases  investigated  by  Rowe 
and  Peaker  (1965) .  In  addition,  this  new  analytical  ap¬ 
proach  does  not  contain  the  physically  inadmissible  over¬ 
lapping  of  the  slip  lines  previously  encountered  in  appli¬ 
cation  of  plasticity  theory  to  such  problems,  nor  does  it 
require  an  assumption  of  stress  discontinuity  within  the 
slip  line  field  as  made  by  other  researchers,  as,  for 
example  Lee  and  Herrington  (1972) .  The  0  values  vary  along 
the  ground  surface,  OA,  indicating  that  the  principal 
stresses  reorient  themselves  as  a  result  of  the  upward 
pushing  of  the  wall.  In  this  procedure,  zone  OCB  disap¬ 
pears,  since  there  is  no  stress  discontinuity  at  0  in  this 


direction. 


CHAPTER  3 


LITERATURE  REVIEW  OF  ANCHOR  PULLOUT  STUDIES 

Baker  and  Konder  (1966)  summarized  the  failure 
hypotheses  for  shallow  anchors  by  the  "friction  cylinder 
method,"  the  "soil  cone  method"  and  Balia's  method.  The 
shapes  of  the  corresponding  failure  surfaces  are  shown  in 
Figure  8.  In  the  friction  cylinder  method,  the  pullout 
capacity,  P,  is  equal  to  the  weight  of  soil,  Wg,  immediately 
above  the  anchor  as  well  as  the  shear  resistance,  S,  along 
the  assumed  vertical  cylindrical  failure  surface.  The  soil 
cone  method  assumes  that  the  failure  surface  takes  the  shape 
of  a  truncated  cone  extending  above  the  anchor  with  a  base 
angle  of  45°  -  <t>/2.  The  pullout  capacity  is  given  as  the 
weight  of  the  soil  within  the  truncated  cone.  The  method 
presented  by  Balia  (1961)  is  based  on  the  shape  of  the  fail¬ 
ure  surface  observed  during  small-scale  model  anchor  tests 
(Figure  8c)  in  sand.  The  pullout  capacity  is  equal  to  the 
weight  of  soil  within  the  assumed  failure  surface  and  the 
side  shear  resistance. 

Prediction  of  the  pullout  capacity  of  a  buried  cir¬ 
cular  plate  is  accomplished  by  the  following  equation: 


*  Assumed  Failure  Surface 


(c) 


Figure  8.  Methods  of  calculating  pullout  capacity  (after 

Baker  and  Konder,  1966).  —  (a)  Friction  cylinder 
method;  { b )  Soil  cone  method;  (c)  Balia's  method. 
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P  =  YDFgA  (10) 

where 

A  =  projected  area  of  the  plate 
D  =  depth  of  the  anchor  plate  below  the  ground 
surface 

F  =  a  pullout  factor  which  is  a  function  of  the 

q 

angle  of  shearing  resistance  of  the  sand,  and 
the  D/B  ratio 
P  =  pullout  capacity 
Y  =  effective  unit  weight 

Modern  research  on  this  subject  began  with  the  paper 
by  Balia.  Baker  and  Konder  (1966)  confirmed  Balia's  major 
findings  regarding  anchor  plates  in  dense  sand;  they  showed 
that  deep  anchors  behaved  differently  from  shallow  anchors. 
Sutherland  (1965)  presented  results  of  pullout  tests  with 
model  plates  up  to  6  in.  in  diameter  in  loose  and  dense 
sand,  as  well  as  with  94  in.  diameter  shafts  in  medium  dense 
to  dense  sands.  He  found  that  the  mode  of  failure  varied 
also  with  sand  density  and  showed  that  Balia's  analytical 
approach  gave  reasonable  results  only  in  sands  of  some 
intermediate  density. 

Vesic  (1971)  showed  that  his  theory  of  expansion  of 
cavities  close  to  the  surface  of  a  semi-infinite  rigid-plastic 
solid,  which  gives  the  ultimate  radial  pressure  needed  to 
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break  out  a  cylinder  or  a  spherical  cavity,  can  be  applied 
to  the  problem  of  anchor  plates.  Vesic  found  that  transi¬ 
tion  to  deep  anchor  behavior  occurred  in  loose  sand  at  D/B 
ratios  of  two  or  three,  and  in  dense  sand  at  D/B  ratios  of 
ten  or  more. 

Healy  (1971)  performed  both  model  and  field  pullout 
tests  in  granular  materials.  He  found  that  the  pullout 
resistance  of  small  anchors  (6  inches)  in  sand  varied  di¬ 
rectly  with  the  depth  provided  that  the  anchors  have  D/B 
ratios  equal  to  or  greater  than  six  in  dense  sand  and  D/B 
ratios  between  two  and  six  in  loose  sands. 

Adams  and  Hayes  (1967)  photographed  pullout  tests  in 
sand  for  D/B  ratios  between  2  and  4.5.  The  failure  surface 
changes  with  relative  density  and  D/B  ratio.  At  shallow 
depths  (D/B  -  2) ,  the  failure  shape  very  closely  approxi¬ 
mates  the  friction  cylinder  method  assumption,  regardless 
of  the  soil  density.  At  greater  depths  (D/B  >  4.5),  the 
failure  zone  was  "local"  for  a  loose  sand  condition,  while 
for  the  dense  sand  tests,  the  failure  zone  reached  the 
ground  surface. 

Kovacs  and  Yokel  (1979)  summarized  the  studies  of 

Bemben,  Kupferman  and  Kalajian  (1971) ,  in  which  the  pullout 

2 

capacity  of  a  55  in  Y  fluke  (marine)  anchor  for  various 
depths  of  embedment  was  presented  as  calculated  by  seven 
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different  theoretical  methods  and  as  measured  in  field 
loading  tests.  There  is  some  agreement  between  predicted 
and  measured  values  for  D/B  ratios  less  than  3.  At  D/B  *  4, 
the  predicted  pullout  capacity  varies  considerably.  As  the 
anchor  is  embedded  further,  only  the  friction  cylinder 
hypothesis  gives  reasonably  conservative  results.  This  was 
also  observed  by  Das  and  Seeley  (1975) .  Based  on  the  wide 
scatter  of  the  pullout  factor,  Fq  versus  D/B  relationship 
as  determined  theoretically  and  in  the  laboratory,  Esquivel- 
Diaz  (1967)  concluded  that  no  satisfactory  theory  is  avail¬ 
able  for  the  determination  of  the  pullout  capacity  of 
circular  disc  anchors. 

The  anchor  pullout  studies  discussed  in  this  chapter 
apply  only  to  dry  sand.  However,  there  are  anchor  pullout 
studies  on  clays  and  marine  soils.  An  extensive  biblio¬ 
graphy  of  anchor  pullout  studies  is  contained  in  Kovacs  and 
Yokel  (1979). 


CHAPTER  4 


SOLUTION  PROCEDURE 

Definition  of  the  Problem 
Figure  9  shows  the  coordinate  axis  used  in  the 
problem.  A  circular  disk  is  embedded  at  some  depth,  +z 
downward,  below  the  surface.  A  force,  P,  acting  along  the 
z-axis  in  the  negative  direction,  is  applied  at  the  center 
of  the  plate.  The  problem  is  axially  symmetric  with  respect 
to  geometry  and  load.  Therefore,  only  half  the  problem 
needs  to  be  solved.  Assumptions  inherent  to  this  study  are: 

1.  Soil  behaves  as  a  rigid-perfectly  plastic  material. 

2.  Soil  strength  properties  are  homogeneous  and  iso¬ 
tropic. 

3.  Seepage  and  inertia  forces  are  equal  to  zero. 

4.  There  are  no  normal  or  shear  stresses  along  the 
ground  surface. 

5.  The  anchor  plate  is  infinitely  rigid  with  respect 
to  the  soil,  so  that  the  soil  fails  before  the 
anchor  plate  deforms. 


The  primary  difficulty  is  determining  the  boundary 
conditions.  It  was  hypothesized  that  as  the  plate  is  moved 
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upward,  a  mass  of  soil  above  the  plate  moves  with  it  as  a 
rigid  body.  The  movement  of  this  rigid  body  of  soil  stres¬ 
ses  the  surrounding  soil,  and  thereby  creates  the  slip  line 
field  at  failure.  The  shape  of  the  rigid  body  is  not  known. 
Two  shapes  were  arbitrarily  selected  and  investigated  separ¬ 
ately  in  this  study.  The  first  was  a  right  circular  cone, 
with  the  apex  at  the  surface  and  the  base  coincident  with  the 
plate  (Figure  9a) .  The  second  was  a  right  circular  cylinder 
extending  vertically  from  the  edge  of  the  plate  to  the  sur¬ 
face  (Figure  9b)  and  having  the  same  diameter  as  the  plate  itself. 


Figure  9.  Assumed  rigid  body  surfaces. 


The  stresses  transferred  from  the  rigid  body  to  the  surround¬ 
ing  soil  are  due  to  passive  pressures,  with  6  <  0  generated 
at  the  soil-rigid  body  interface.  Therefore,  the  analytical 
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approach  of  Nowatzki  and  Karafiath  (1981)  was  extended  to 
the  problem  at  hand,  with  the  rigid  body  surf  act1  acting 
similar  to  a  retaining  wall  boundary. 

Procedure 

Nowatzki  and  Karafiath' s  computer  program  was  al¬ 
tered,  and  the  boundary  conditions  outlined  in  the  previous 
section  were  incorporated.  The  program  for  each  rigid  body 
is  included  in  Appendix  A.  Figure  10  shows  the  numbering  of 
nodal  points  and  the  layout  of  the  problem  for  computer 
solution.  The  solution  procedure  is  demonstrated  schematic¬ 
ally  in  Figure  11.  Equations  9  were  used  to  determine  x  and 
0  values  along  the  ground  surface,  OA,  and  Equations  8  were 
used  to  compute  x,  z,  o,  and  0  values  of  the  remaining  nodal 
points.  Figure  10  also  shows  the  plot  of  6  values  along  the 
ground  surface,  with  negative  0  values  shown  upward. 

Although  z  =  0  and  a  =  0  along  the  ground  surface, 
OA,  an  explanation  is  required  of  how  the  parameters  were 
obtained  along  the  rigid  body  surface,  OD.  After  the  soil 
parameters  and  the  geometry  of  the  problem  are  input,  the 
program  divides  the  boundary  OD  into  (n  -  1)  equal  segments 
(refer  to  Figure  10) ,  where  n  =  some  convenient  number  of 
nodal  points.  For  the  cylindrical  rigid  body  surface,  x  =  0 
and  Az  =  depth  of  the  plate/ (n  -  1) .  For  the  cone  surface. 
Ax  =  plate  radius/(n  -  1)  and  Az  =  depth/(n-l) . 
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RigiS  Body  Surface; 
Ground  Surface;  z  « 


x,  z,  □,  and  6  known 
a  »  0 ;  x  and  9  unknown 


Figure  11.  Solution  procedure. 


Along  both  surfaces,  6  was  assumed  to  remain  constant 
and  the  surface  load  varied  linearly  from  zero  at  the  sur¬ 
face,  point  0,  to  the  maximum  value  at  the  base,  point  D. 

The  stress  at  this  point  was  designated  q^.  The  magnitude 
of  6  and  q^  was  estimated  and  input  for  an  initial  trial. 
Then,  6  and  a  at  point  D  was  determined  using  the  following 
relationships ; 
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( 1/2) arctan 


sin  ( 


sin6 

sin$ 
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For  the  rigid  body  cone,  the  inclination  of  the  rigid  body 
surface  (i.e.,  angle  3  of  Figure  9a)  had  to  be  taken  into 
account  in  determining  9. 

The  above  procedure  required  assuming  various  combi¬ 
nations  of  6  and  values,  because  many  of  these  assumed 
combinations  of  values  did  not  create  a  proper  slip  line 
field  (i.e. ,  overlaps  occurred  within  the  slip  line  field) . 

A  solution  was  assumed  to  be  acceptable  when ,  for  a  given  qb ,  6 
was  varied  until  ocerlapping  in  the  slip  line  field  was  elimina¬ 
ted.  The  acceptable  combination  of  6  and  q^  that  predicted  the 
greatest  pullout  capacity  was  assumed  to  be  the  correct  one. 

Figures  12  through  15  demonstrate  the  variations  in 
the  slip  line  field  with  different  assumed  6  and  qfc  along 
the  rigid  body  surfaces.  These  figures  also  show  the  vari¬ 
ations  in  0  along  the  ground  surface.  In  each  of  these 
figures  the  middle  plot  was  assumed  to  be  the  correct  slip 
line  field.  To  produce  Figures  12  and  13,  6  value  remained 
constant  and  the  surface  loading  was  progressively  increased 
from  (a)  to  (c) .  Overlaps  within  the  slip  line  field  are 
circled.  For  Figures  14  and  15,  the  surface  loading  re¬ 
mained  constant  and  6  values  were  varied,  becoming  increas¬ 
ingly  more  negative  from  (a)  to  (c) .  Note  that  Figures  14a 
and  15a  are  also  acceptable  slip  line  fields.  However, 
neither  is  considered  to  be  the  correct  one,  in  accordance 
with  assumptions  stated  in  previous  paragraphs. 
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The  process  of  assuming  various  combinations  of  6 
and  and  evaluating  the  results  with  regard  to  permissible 
solutions  was  greatly  facilitated  by  using  the  facilities 
of  the  University  of  Arizona  Interactive  Graphics  Engineer¬ 
ing  Laboratory.  Here  the  x,  z,  o,  and  0  values  for  the 
j  =  n  nodal  points  were  viewed  on  the  screen  as  they  were 
calculated,  i.e.,  in  real  time.  Generally,  if  there  was  an 
overlap  within  the  slip  line  field,  it  occurred  on  the  j  =  n 
characteristic  line.  Plotting  was  then  done  as  a  check  on 
what  otherwise  seemed  an  acceptable  slip  line  field. 
Therefore,  this  iterative  procedure  was  easily  accomplished 
internal  to  the  execution  of  the  program.  This  procedure 
saved  time,  eliminated  redundant  calculations,  and  decreased 
the  possibility  for  human  error. 

The  CALPLOT  program  was  used  to  plot  the  mesh  of  x 
and  z  values  calculated  and  the  0  values  along  the  ground 
surface.  The  CALPLOT  program  consists  of  FORTRAN  callable 
subroutines,  each  of  which  works  through  the  computer  hard¬ 
ware  system  to  drive  the  Calcomp  602  plotter.  A  sample 
printout  of  slip  line  field  values  is  included  in  Appendix  B. 

Once  the  correct  slip  line  field  was  determined,  the 
pullout  capacity,  P,  and  pullout  factor,  F  ,  were  calculated 

4 

as  follows.  The  z-direction  component  of  q^  was  designated 
as  o  .  Integration  of  the  triangular  stress  distribution 
over  the  rigid  body  surface  yielded  the  following  expressions 
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For  the  cone : 

S  =  (1/2)  o2  tt  R  /r2  +  D2 


(12a) 


For  the  cylinder: 

S  =  a  ttR  D 
z 


(12b) 


Adding  the  shear  resistance,  S,  to  the  weight  of  the  soil 
within  the  rigid  body  yielded  the  pullout  capacity  as 
follows : 

For  the  cone : 

P  =  S  +  (1/3)yttR2D  (12c) 

For  the  cylinder: 

P  =  S  +  yitR2D  ( 12d) 


Rewriting  Equation  10  yielded  the  following  expression  for 
pullout  factor: 


where 


Fq 


2 

YTtR  D 


( 12e) 


R  =  plate  radius  (B/2) 
D  =  plate  depth 
Y  =  soil  unit  weight 


The  pullout  factor  is  used  in  many  studies  to  com¬ 
pare  pullout  capacity  results  with  other  reports.  In  this 
way,  different  sized  anchors  can  be  equated  if  their  relative 
depth  (D/B)  and  the  soil  parameters  are  the  same. 


CHAPTER  5 


PRESENTATION  AND  DISCUSSION  OF  RESULTS 

A  summary  of  the  numerical  results  for  this  study 

is  given  in  Table  1.  The  relationship  between  pullout 

factor  (F  )  and  relative  depth  (D/B)  determined  in  this 
Si 

study  are  plotted,  along  with  other  theoretical  and  experi¬ 
mental  results  in  Figures  16a  and  16b  for  loose  and  dense 
sand,  respectively.  In  this  study,  loose  sand  was  consid¬ 
ered  to  have  a  unit  weight  of  100  pcf  and  an  angle  of  in¬ 
ternal  friction  of  31°;  dense  sand  was  assumed  to  have  y  = 

112  pcf  and  <J>  =  42°. 

The  theoretical  results  of  this  study  are  lower  than 
experimental  findings,  but  are  much  closer  at  shallow  depths, 
with  the  differences  becoming  progressively  greater  with 
depth.  The  rigid  body  cylinder  assumption  predicts  higher 
pullout  capacities  than  the  cone.  A  numerical  comparison 
of  the  theoretical  F^  values  with  experimental  ones  is  shown 
in  Table  2.  For  D/B  ratios  down  to  four,  the  cylinder  F^ 
values  in  loose  sand  are  an  average  of  81  percent  of  Suther¬ 
land's  results.  At  D/B  =  1,  theoretical  values  are  in 
agreement  with  Vesic's  Theory,  but  rapidly  diverge  with 

depth.  For  dense  sand,  the  cylinder  F  values  are  an  average 
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Rigid  Body 

D/B 

Ratio 

Angle  of 
Internal 
Friction  (<(>) 

Unit 
Weight 
(pcf ) 

Pullout 

Capacity’ 

(lb) 

* 

F 

q 

Cone 

1 

31 

100 

1.72 

1.  39 

II 

2 

It 

It 

4.86 

1.98 

If 

3 

II 

II 

9.56 

2.59 

II 

4 

It 

It 

15.57 

3.17 

ft 

10 

II 

II 

Cylinder 

1 

It 

II 

2.63 

2.14 

It 

2 

11 

ft 

8.07 

3.29 

II 

3 

It 

II 

16.33 

4.44 

II 

4 

It 

It 

27.36 

5.57 

M 

10 

II 

It 

152.94 

12.46 

Cone 

1 

42 

112 

2.35 

1.70 

II 

2 

II 

It 

7.05 

2.56 

M 

3 

II 

It 

14.13 

3.43 

It 

4 

II 

It 

23.57 

4.28 

II 

10 

It 

It 

129.33 

9.41 

Cylinder 

1 

If 

It 

3.73 

2.72 

II 

2 

It 

It 

12.19 

4.44 

It 

3 

II 

II 

25.36 

6.14 

ft 

4 

II 

It 

43.28 

7.87 

•1 

10 

II 

II 

249.84 

18.18 

*  For  3"  diameter  plate 


Pullout  Factor,  Fa  Pullout  Factor 


of  58  percent  of  Baker  and  Konder's  and  40  percent  of 
Sutherland's  results. 

Figure  17  shows  that  theoretical  pullout  factors 
increase  linearly  with  relative  depth  increase.  In  contrast, 
experimental  values  appear  to  increase  exponentially  with 
increasing  relative  depth. 

The  causes  of  the  disagreement  with  observed  data 
were  not  determined.  However,  the  assumption  of  linearly 
varying  a  and/or  constant  0/6  along  the  surface  of  the  rigid 
body,  could  be  refined.  Varying  these  parameters  along  the 
rigid  body  surface  might  produce  an  acceptable  slip  line 
field  at  a  more  inclined  6,  thus  mobilizing  more  shear 
strength  of  the  soil  and  increasing  the  predicted  pullout 
capacity.  Perhaps  another  rigid  body  shape  is  correct, 
possibly  a  cone  truncated  at  the  ground  surface,  or  an 
elliptical  or  parabolic  rigid  body  surface.  Investigation 
of  these  shapes  was  beyond  the  scope  of  this  study. 

It  should  be  noted  that  the  plane  strain  equations  were 
used  in  this  analysis,  and  not  the  axially  symmetric  ones .  The 
latter  equations  are  presented  in  Karafiath  and  Nowatzki  (1978)  . 

Attempts  were  made  to  use  the  recurrence  relation¬ 
ships  for  the  axially  symmetric  case;  however  the  solutions 
did  not  converge.  It  is  felt  that  the  boundary  conditions 
imposed  by  the  anchor  pullout  problem  made  the  recurrence 
relationships  ill-conditioned  for  the  numerical  method 
followed  in  the  plane  strain  solution  procedure. 
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Comparison  with  Coulomb 
Trial  Wedge  Method 

Since  not  all  engineers  have  ready  access  to  a 
computer,  it  was  desired  to  find  an  easy  hand  method  to 
approximate  the  theoretical  analysis.  It  was  found  that 
the  results  of  the  Coulomb  graphical  trial  wedge  method  (see 
Figure  18)  correlated  well  with  the  computer  results.  For 
6  =  0.750  as  obtained  from  the  computer  solution,  the  most 
critical  failure  wedge  with  the  rigid  body  surface  was  found 
to  be  approximately  30°  for  both  the  loose  and  dense  soil. 

The  Z-directional  component  of  the  measured  force  was 
integrated  over  the  circumference  of  the  rigid  body  cylinder 
and  then  added  to  the  weight  of  the  soil  within  the  rigid 
body  cylinder  to  determine  the  pullout  capacity,  P.  The 
pullout  factor,  F  ,  was  then  calculated  using  Equation  12e. 
There  is  a  very  close  correlation  between  the  results  of 
the  computer  analysis  and  those  of  the  trial  wedge  method. 

The  comparison  is  shown  in  Table  3. 

Thus  it  seems  that  the  graphical  solution  could  be 
used  just  as  well  as  the  sophisticated  computer  solution 
to  obtain  an  estimate  of  Pullout  Factor  F  for  the  conditions 
and  assumptions  made  in  this  study.  Obviously,  both  solu¬ 
tions  fall  short  of  the  test  data,  indicating  that  the 
plane  strain  assumption  may  be  unrealistic  for  this  case. 


Figure  18.  Trial  wedge  method.  —  (a)  Failure  wedge  and 
acting  forces;  (b)  force  polygon  for  compu¬ 
tation  of  P  . 


Table  3. 

Correlation  between  computer 
wedge  method  Fg  values. 

analysis  and  trial 

D/B 

<p  Computer 

Trial  Wedge 

1 

31° 

2.14 

2.17 

2 

31° 

3.29 

3.33 

3 

31° 

4.44 

4.50 

4 

31° 

5.57 

5.67 

1 

42° 

2.72 

2.75 

2 

42° 

4.44 

4.50 

3 

42° 

6.14 

6.24 

4 

42° 

7.87 

7.99 

CHAPTER  6 


CONCLUSIONS 

This  has  been  a  theoretical  anchor  pullout  study, 
using  the  plasticity  theory  and  hypothesis  that  a  rigid  body 
of  soil  moves  with  the  anchor  plate  as  it  is  moved  toward 
the  surface.  Two  such  rigid  body  shapes  were  investigated: 
the  right  circular  cylinder  and  right  circular  cone.  The 
theoretical  pullout  capacities  were  found  to  be  lower  than 
experimental  results,  with  the  cylinder  rigid  body  provid¬ 
ing  closer  agreements  than  the  cone.  At  shallow  depths  in 
loose  sand,  predicted  pullout  capacities  are  approximately 
81  percent  of  experimentally  determined  ones.  In  dense 
sand,  predicted  values  are  about  50  percent  of  experimental 
ones. 

Assumptions  of  this  report  might  be  changed  to  pro¬ 
vide  better  results.  Varying  6  and/or  applying  a  non- 
triangular  stress  distribution  along  the  surface  may  produce 
acceptable  slip  line  fields  at  a  more  inclined  6,  thus  pre¬ 
dicting  a  greater  pullout  capacity.  Additionally,  there  is 
also  the  possibility  that  neither  of  the  rigid  body  surfaces 
investigated  are  the  correct  ones. 
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A  conventional,  graphical  solution  procedure  was 
also  used  to  determine  the  pullout  capacity.  The  results 
of  that  solution  correlate  well  with  the  computer  results. 
This  suggests  that  small  engineering  offices,  without 
access  to  a  computer,  may  be  able  to  utilize  simpler, 
approximate  solutions  for  anchor  pullout  problems  provided 
the  conditions  assumed  in  this  study  pertain. 

Recommendations  for  Further  Study 

1.  Refine  the  mathematical  model  to  account  for  vari¬ 
able  6  and/or  variable  stress  distribution  along  the 
rigid  body  surface. 

2.  Investigate  other  possible  rigid  body  surfaces, 
such  as  a  truncated  cone,  ellipse,  or  parabola. 

3.  Perform  full-scale  and  model  tests. 

4.  Further  study  of  the  axially  symmetric  computer 


solution. 


APPENDIX  A 

COMPUTER  PROGRAMS 
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CONE  RIGID  BODY 

D I  HENS  I  ON  X  (SI.  U).  Z(21,  U),SIC(21.  U).THE(21.  11) 

OPEN  12. "ARRAYS.  OUT",  ATT-"SOM" 

ACCEPT  "INPUT  PHI, GAMMA-  ".PHI, GAMMA 

TYPE  "  PHI-", PHI."  GAMMA-".  GAMMA 

ACCEPT  "INPUT  DEPTH, LENGTH  “, DEPTH. HLENG 

TYPE  "DEPTH-", DEPTH,  "  LENGTH-".  HLENG 

ACCEPT  "INPUT  N  ",  N 

TYPE  “N— " , N 

GO-GAMMA 

ARC-. 0174533 

PI-3  14159 

F1-ARC*PHI 

TF«TAN(F1 > 

U1-PI/4-F1/2 
DO  58  J-l. N 
DO  58  I  — 1 .  21 
X(I, J)-0 
Z(  I.  J)=  0 
SIGU.  U)=0 
THE(I.  J)-0 
NT-N 

ACCEPT  "INPUT  DELTA,  GBOT  ".DELTA. QBOT 
TYPE  "DELTA-".  DELTA,  "  QBOT-",  OBOT 
DEL=ARC*DELT  A 
BETA-ATAN ( HLENG /DEPTH ) 

DTHE-.  5* <  ATAN<SIN<SIN<DEL  > /SIN (FI > )/ 

SORT  < l-( SIN (DEL) /SIN (FI ) >**2) )+DEL)-BETA 
S I BOT— GBOT/ <  COS  <  DEL ) +SGRT  <  COS ( DEL ) **2-C0S (FI) **2  > ) 

DO  80  J-l.N 
I-N+l-J 

F- ( FLOAT ( U ) - 1 ) / <  N— 1 ) 

X(I, J)-F*HLENG 
ZU.  J)-F*DEPTH 
SIGU.  J)-F*SIBOT 
THE(I. J ) -DTHE 

IF (ABS(THE(N, 1 ) ) . GT. U1 >  NT-N-1 
TYPE  "THE, U1-".  THE(N,  1 ),  U1 
DO  90  J-2.  N 

DO  90  I-(N+2-J), (NT+J-l) 

IF< I. EG.  (NT+J-l ) )  GO  TO  10 
L— 0 

TH1-THE( I.  J-l )+Ul 
TH2-THEU-1,  J)-U  1 

511- SIGU,  J-l  ) 

512- SIG ( 1-1 , J> 

V7=2*SICU-1,  J)*SIG<  I.  J-l  ) 

U9=SIG<  I  —  1 ,  J)+SIGU,  J-l  ) 

V8- ( THE ( I , J-l )-THE( 1-1.  J) )*TF 
V9=2*5IG< 1-1 ,  J)*SIG(I,  J-l >*VB 

U6-2*TF*(S1G(I,  J-l )*THE( I.  J-l >+SlG< 1-1, J)*THE( 1-1,  J) ) 
V1-TAN(TH1 ) 

V2-TAN( TH2 ) 

XI-V1+X< I. J-l ) 

X J-V2*X  < I  —  1 ,  J) 

VI 2=1 / ( V1-V2 ) 

X(I,  J>— V12*(Z<  1-1.  J)-Z(I,  J-D+XI— XJ> 

Z(I,J>*ZU-1,J)  +  (XU,J)-X(I-1,  J)  )*V2 
AA--TF* ( X ( I ,  J)-X(I-1,  J) >  +  (Z(I.  J)-Z( 1-1,  J) ) 

EB-TF* ( X  < I ,  JI-XII,  J-l ) )  +  (Z(I.  J)-Z( 1,  J-l ) ) 

U5-SIGU,  J-l  ) -SIGU- 1.  J) 

I F ( U9  EG  0)  GO  TO  15 

SIGU,  J )  -  ( V7+V9+G0*  ( S 1 1  *AA+SI2*BB  )  ) /U9 
THEU.  J )  -  <  U5+U6+G0*  ( BB-AA  )  )/<2*TF#U9) 
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GO  TO  16 

15  SIG< I ,  J)“GO* ( AA+BE ) / ( 2+TF* (THE  < 1-1 • J)-THE( I. J-l) >  > 

THE ( 1 . J)«(BB*( 1+THE( 1-1.  J) )-AA*( 1-THE(I,  vl-1 > > )/(AA+BB> 

16  IF(L. EQ.  1 )  GO  TO  20 

TH1  =  5*  < THE ( I . J-l )+THE< I. J) )+Ul 
TH2«.  5* ( THE ( I  —  1 . JJ+THEC I.  J> )-UI 
SI1*.  5*(SIG< I.  J-l )+SIG< 1. J) ) 

SI2=.  5*(SIG(I-1,  J)+SIG(I, J>> 

U6«2*TF*(SI1*TH£! I,  J-l )+SI2*THE ( 1-1 , J) ) 

V7*S11*SIG< I— 1 . J)+SI2*S1G(I, J-l > 

V9*2+SI1*SI2*VS 
U9«SI 1+SI2 
L=  1 

GO  TO  5 

10  ZCI,J>=0 
SIG  < I > J>«0 

TH3=T  AN ( THE ( I - 1 . J )  -U 1 ) 

KK=0 

11  X( I, J)«X< I  — 1 . J)-Z( 1-1.  J)/TH3 

AA=-TF* ( X ( 1  /  J)-X(I-1»  J) >  + (Z<1. J)-Z( 1-1. J) ) 

U3=SIG ( I , J)-S1G< 1-1. J>-CO*AA 

THE < I .  J)=U3/(2*SIG< 1-1.  J)*TF>+ THE(I-1.  J> 

TH3=TAN( .  5*(THE(  I.  J)+THE<1-1,  J>  )-Ul  > 

IF (KK.  EG.  6)  GO  TO  20 

KK-KK+1 

GO  TO  11 

20  IF<J.  EG.  N>  TYPE  “I.J,  X,  Z-  ”,  I ,  J,  X  <  1 ,  J  > .  Z  <  I ,  J) 

90  CONTINUE 

S I ZmQBOT*S 1 N < BETA-DEL ) 

S“.  5*SX  Z*P I *HLENG*SGRT (HLENG**2+DEPTH**2  > 

P-S+GO*P I *HLENG*+2*DEPTH/3 
FG=P / < GO*DEPTH*P I *HLENG**2 > 

TYPE  "P»",P,"  FG«".FG 

ACCEPT  "INPUT  1  FOR  CHANGE  OF  DATA, ELSE  0  ",IJK 
IF< IJK. EG.  1 )  GO  TO  59 

ACCEPT  "INPUT  1  FOR  PRINT  OF  DATA, ELSE  0  ",KJI 
IF (KJI .  EG  0)  GO  TO  75 
WRITE  <12.24) 

24  FORMAT  </39X, "CONE  RIGID  BODY" / "PHI , CO, DEPTH, HLENG. QBOT, " , 
1  "P.FQ, DELTA") 

WRITE  (12.25)  PHI. GO,  DEPTH, HLENG,  QBOT,  P,  FG.  DELTA 

25  FORMAT  (10F10  4) 

WRITE  (12,30) 

30  FORMAT  < //40X, "X-COORDINATES" // > 

WRITE  (12-  35) <  <X< I, J), J-l, N),  1*1,  (NT+10) ) 

35  FORMAT  (11F11.4) 

WRITE  (12,40) 

40  FORMAT  C//40X. " Z-COORD INATES" /V  > 

WRITE  ( 12,  35) ( <Z< I,  J),  J*l,  N),  I«I,  (NT+10) ) 

WRITE  (12,50) 

50  FORMAT  < //40X, "SIGMA  VALUES" /f ) 

WRITE  < 12, 35) <  <SIG< I. J), J»l,  N).  1*1.  (NT+10) ) 

WRITE  (12,60) 

60  FORMAT  ( //40X, "THETA  VALUES"//) 

WRITE  (12, 35) ( <  THE  < I,  J).  J*l,  N),  1-1.  (NT+10) ) 

75  ACCEPT  "INPUT  1  FOR  PLOT,  ELSE  0  ",M 
ACCEPT  "INPUT  SCALE  ",SC 

IF(M  EQ  1)  CALL  PLTSL(X,  Z, THE,  N, NT, HLENG,  DEPTH,  SC) 

CLOSE  12 

STOP 

END 
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SUBROUTINE  PLTSL(X,  Z,  THE,  N,  NT,  HLENO,  DEPTH,  80 
DIMENSION  X (21,  11),  2(21,  ll),THE(2l. 11) 

CALL  1 14602  ( 1,  "PASS.  PLT*  ) 

DO  6  U-l.N 

DO  6  1-1. (2*N-1) 

X<I,  J)-SC*X(I,  J> 

6  Z(I,  J>  —  SC*Z(I,  J> 

CALL  PLT602  <7.  9,0.  ,  3) 

CALL  PLT602  (O.  ,0.  ,2> 

CALL  PLT602  (  0.  , (-SC*DEPTH) . 2 > 

CALL  PLT602  ( <SC*HLENG> ,  (-8C*DEPTH> ,  2> 

CALL  PLT602  (  0.  ,0.  .2) 

DO  9  J-2.  N 
K-N+l-J 

CALL  PLT602  ( X(K,  J) ,  ( Z (K, J) ) ,  3) 

DO  9  I«(K+I  ),  (NTt-J-1 > 

3  CALL  PLT602  (X(I,  J),  (Z(X,  J>),2> 

DO  3  1-2. (NT+N-2) 

K-N+l-I 

IFU.OT.  N)  K-I-NT-H 

CALL  PLT602  ( X ( I,  K> .  <  ZU,  K>  > ,  3) 

DO  3  U-(K+1>.N 

3  CALL  PLT602  (X(I,  J),  (Z(X,U>>,2> 

CALL  PLT602  (7.  9, .  29,  3) 

CALL  PLT602  (0.  ,  .  29.  2) 

CALL  PLT602  ( X(N. 1 ) ,  (-. 9*THE(N,  1 >+.  29>.  2) 

DO  2  I-N,  (NT+N— 1 ) 

U-I-NT+1 

2  CALL  PLT602  (X(I,  J),  (-.  9*THE( I.  J>+.  29),  2) 
CALL  PIN602 
RETURN 

END 
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CYLINDER  RIGID  BODY 

DIMENSION  X<21.  11 >. Z<21.  11 >, SXG(21.  11 ).  THE <21.  11 ) 

OPEN  12.  "ARRAYS.  OUT".  ATT-"SO M" 

ACCEPT  "INPUT  PHI.  GAMMA-  ",  PHI.  GAMMA 

TYPE  "  PHI-", PHI,"  GAMMA-",  GAMMA 

ACCEPT  "INPUT  DEPTH, LENGTH  ",  DEPTH,  HLENG 

TYPE  "DEPTH-". DEPTH, "  LENGTH-" , HLENG 

ACCEPT  "INPUT  N  ",N 

TYPE  "N=",N 

GO-GAMMA 

ARC-.  0174533 

PI-3. 14159 

F1=ARC*PHI 

TF-TAN(F1 ) 

U1-PI/4-F1/2 
59  DO  58  U-l-N 
DO  58  1=1,21 
X ( I , J ) =0 
2  < I , J)=0 
SIG< I, J)=0 
58  THE  < I ,  J ) =0 

NT-N 

ACCEPT  "INPUT  DELTA,  GBOT  ",  DELTA,  GBOT 
TYPE  "DELTA-", DELTA. "  GBOT-", QBOT 
DEL=ARC*DELTA 

DTHE-.  5*  <  AT AN (SIN(SIN( DEL )/SIN(Fl) )/ 

1  SQPT ( 1- (SIN  <  DEL > /SIN(F1 ) >**2> >+DEL) 

S IBOT-GPOT/ <  COS  < DEL  >+SGRT  <  COS ( DEL ) *#2-C0S  <  FI >**2 ) > 

do  eo  J-l. N 

I-N+l-J 

F-<FL0AT!U)-1)/<N-1> 

X  < I, U> -HLENG 
Z< I, U ) —F* DEPTH 
SIG ( I ,  J)«F*SIBOT 
BO  THE< I, U3-DTHE 

IF  <  ABS ( THE  < N,  1  >  ) .  GT.  U1 )  NT-N-1 
TYPE  "THE, U1-  ", THE(N,  1 ) ,  Ul 
DO  90  J-2,  N 

DO  90  I-CN+2-U),  <NT+U-1 > 

IF< I. EQ.  (NT+U-l ) )  GO  TO  10 
L— 0 

TH 1 —THE  < I , U- 1 >  +U 1 
TH2=THE< 1-1 , U3-U1 
SI 1=SIG< I, U-l ) 

SI2-SIG ( I — 1 . U> 

V7=2*SIG(I-1,  U)«SIGU.  U-l) 

U9=S1G< I  — 1 ,  U)+SIG< I, U-l ) 

VS—  t THE  < I , U-I ) -THE ( I - 1 , U) )*TF 
V9*2*SIG ( I- 1 • U)*SIG< I,  U-l >*VB 

U6=2*TF*CSIG 'I,  U-l ) *THE ( I ,  U-l J+SIGCI-l,  U)*T HEC 1-1, U) ) 
5  VI— TAN  <  TH1 : 

V2— TAN  f  TH2  > 

X I -VI *X ( I , U-l > 

XU— V2* X  < 1-3  ,  U) 

V 1 2- 1 / <  V 1 - V2 ) 

X  <  I ,  UJ-V12-*  (  Z  <  I-l .  U  )-Z  ( I ,  U-l J+XI-XU) 

Z(I, U)-Z< 1-1 . U)  +  (X< 1. U)-X< I- 1,  U) > *V2 
AA— TF#  <  X  <  J  ,  U)-X(  I- 1,  U)  )  +  <Z<  1,  U)-Z(I-1,  U)  ) 

BB-TF* ( X ( 1 ,  U)-X< I,  U-l ) )  +  <Z< I. U>-Z(l, U-l ) ) 

U5-SIG  < I , U-I )-SIG< I - 1 . U) 

IF(U9  ECi.  0)  GO  TO  IS 

SI G ( 1 , U ) - <  V7+V9+G0*  <  SIJ  *AA+S1 2*EE ) ) /U9 
THE  f I ,  U ) - <  U5-*-U6>+GO*  <  BB-AA ) > / C  2*TF*U9 ) 

GO  TO  16 


53 


15  S1GU,  J  >=GCI*  ( AA+BE  1  /  (2+TF* (THE (  3  —  1 »  J)-THE(I.  J-l > )  1 
THE(I,  J)=(BB*(1+THE(I-1,  J> >-AA* ( 1-THE ( I ,  J-l ) > ) / 1 AA+BB > 

16  1F(L. EQ.  11  GO  TO  SO 

TH1*.  5*  ( THE ( I .  J-l  >+THE( I.  J>  1+U1 
TH2*.  5* (THE (1-1.  J)*THE(l,  J) >-Ul 
SI l*. 5*(SIG(1.  J-1)+SIG(I. J) 1 
SIS*.  5*(SIG(I-i. J)+£IG< I. J) ) 

U6=2*TF* (SI 1*TH£( I .  J-l 1 +SI2*THE ( 1-1 . J) 1 
V7*5I 1*SIG( 1-1. J)+S12*SIC( I, J-l> 

V9=S*SI1*S12*V8 

U9*SI1+SI2 

L*1 

GO  TO  5 

10  Z(I,J)=0 
SIG( I. J!*0 

TH3*T  AN ( THE ( I - 1 , Jl-Ul ) 

KK=0 

11  XtJ, J)*X(I-1. J)-Z(I-1. J1/TH3 

AA«-TF  *  (  X  (  I .  J)-X(I-t,J)>+(Z(I,J>-Z(I-l,J>> 

U3»SIG( 1 . J)-SIG( I-l.  J)-GO*AA 

THE (I.  J)*U3/ <2*SIG< t-1, J)*TF)+THE(I-1.  J) 

TH3*TAN(.  5* ( THE ( I . J)+THE(X-1, J>  >-Ul ) 

I F ( K.K .  EQ  6)  GO  TO  20 

KK-KK+1 

GO  TO  11 

20  IF(J.EQ.N)  TYPE  "I,  J,  X.  Z*  ",I,J,  X(  l.Jl.ZU.J) 

90  CONTINUE 

SIZ*-QE0T*S1N(DEL) 

P=S I Z*P I *HLENG*DEP TH+GO*P  t  *HLENG**2#DEPTH 
FQ*P / ( GO*DEP  TH*P 1 *HLENG**2 ) 

TYPE  "P*",P,"  FQ*'\FQ 

ACCEPT  "INPUT  1  FOR  CHANGE  OF  DATA.  ELSE  O  ".IJK 
IF< I JK  EQ.  1 )  GO  TO  59 

ACCEPT  "INPUT  1  FOR  PRINT  OF  DATA,  ELSE  0  ",KJI 
tF<KJI.  EQ  0)  GO  TO  75 
WRITE  (12.24) 

24  FORMAT  ( /37X. "CYLINDER  RIGID  BODY" .  /"'PHI , GO, DEPTH. HLENG. QBOT,  " 

1  , "P, FQ. DELTA" ) 

WRITE  (12,25)  PHI. CO, DEPTH.  HLENG,  QBOT,  P, FQ, DELTA 

25  FORMAT  (10F10  4) 

WRITE  (12.30) 

30  FORMAT  (//40X. "X-CQQRDINATES'V/ 1 

WRITE  (12. 35) ( (XII. J>. J*1,N>, 1*1.  (NT+10) > 

35  FORMAT  ( 1 IFl 1 . 4 ) 

WRITE  (12,40) 

40  FORMAT  ( //40X, "Z-CDORDINATES" / / 1 

WRITE  <12.  35) < (Z(I,  J),  J*l,  N), 1*1,  (NT* 10) > 

WRITE  (12.501 

50  FORMAT  <//40X, "SIGMA  VALUES"//) 

WRITE  ( 12, 35) ( (S1C( I, J), J-l, N). 1*1- (NT+10) 1 
WRITE  (12. 60) 

60  FORMAT  (//40X. "THETA  VALUES"// 1 

WRITE  (12, 35' ( ( THE ( I , J) , J*1 , N) , I-l , (NT+10) ) 

75  ACCEPT  "INPUT  1  FOR  PLOT,  ELSE  0  ",M 
ACCEPT  "INPUT  SCALE  ",  SC 

IPdl.EO  1)  CALL  PLTSLFU,  Z.  THE,  N,  NT,  HLENG,  DEPTH,  SC) 

CLOSE  12 

STOP 

END 
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SUBROUTINE  PLTSLF ( X.  Z.  THE.  N,  NT.  HLENC,  DEPTH.  SC ) 
DIMENSION  X (21. 11).  Z(21.  11), THE (21. 11) 

CALL  IN602  ( 1. "PASS.  PLT" ) 

DO  6  J-l.N 

DO  6  I-l.  <2«N-1) 

X(I. J)-SC*X<I,  J> 

Z(I. J>  — SC*Z(I.  J> 

CALL  PLT602  (7.  9.0.  ,3) 

CALL  PLT602  (0.  .0.  .2) 

CALL  PLT602  (  0.  ,  <-SC*D£PTH> ,  2) 

CALL  PLT602  <(BC*HLENC>,  <-SC*DEPTH>,  2) 

CALL  PLT602  ( <SC*HLENQ> ,  0.  .  2) 

DO  9  J-2.N 
K-N+l-J 

CALL  PLT&02  <X(K,  J).  <Z(K.  J>).3) 

DO  9  I-(K+1). (NT+J-l) 

CALL  PLT602  < X< I. J) .  < Z < I.  J) >.  2) 

DO  3  1-2,  (NT+N-2) 

K-N+l-I 

IFU.OT.  N)  K-I-NT+1 

CALL  PLT602  ( X ( I. K ) ,  ( Z( I, K ) ) .  3) 

DO  3  J-<Ktl),N 

CALL  PLT602  (X ( I. J) .  < Z< I, J) >,  2) 

CALL  PLT602  <7.  9.  .  29.  3) 

CALL  PLT602  ( < X (N.  1 ) ) .  .  29. 2) 

CALL  PLT602  <X(N.  1 ).  (-.  9*THE(N. 1 >♦.  29).  2) 

DO  2  I-N.  (NT+N-l) 

J-I-NT+1 

CALL  PLT602  (X< I. J) .  9*THECI.  UK.  29) .  2) 

CALL  FIN602 

RETURN 

END 
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